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Abstract 

We compare the distortion effects of three popular final-state hadronization schemes. We show 
how hadronization modifies the initial-state gluon correlations in high energy p-i-A collisions. 
The three models considered are (1) LPH: local parton-hadron duality, (2) CPR; collinear parton- 
hadron resonance independent fragmentation, and (3) LUND: color string hadronization. The 
strong initial-state azimuthal asymmetries are generated using the GLVB model for non-abelian 
gluon bremsstrahlung, assuming a saturation scale Qsat = 2 GeV. Long-range elliptic and triangu¬ 
lar harmonics for the final hadron pairs are compared based on the three hadronization schemes. 
Our analysis shows that the process of hadronization causes major distortions of the partonic az¬ 
imuthal harmonics for transverse momenta at least up to pr - GeV. In particular, they appear 
to be greatly reduced for < 1 2 GeV. 

Keywords: 

Heavy Ion Collisions, QCD Radiation, Monte Carlo, Hadronization 
PACS: 24.85.H-p, 25.75.-q, 21.60.Ka 


1. Introduction and motivation 


Multi-particle, long-range in pseudo-rapidity azimuthal correlations are widely studied in 
high-energy nuclear collisions at RHIC Auh-Au |[T]|2l and LHC Pb-i-Pb |[3}j5l- In particular, they 
are considered as a signature for the “perfect fluid” behavior of the strongly coupled Quark- 
Gluon-Plasma (sQGP) produced in such reactions. 

The recent discovery of long-range pH- A azimuthal harmonics, see Pig. 1 1(a) with magnitudes 
and pt dependence comparable to the ones found in A-i-A collisions M, and the near energy 
independence of these A-hA moments observed in the Beam Energy Scan (BES) at RHIC Col 
together with LHC have challenged the uniqueness of the sQGP interpretation of v„ in Ah-A. Prior 
to the recent p-i-A and BES data, it was assumed that in smaller transverse size p-i-A systems 
or lower energy Ah-A reactions, the perfect fluid sQGP “core” would gradually change into a 
highly dissipative hadron resonance gas “corona” and lead to a substantially smaller magnitude 
of the azimuthal harmonic moments. The observed near independence of Vnipr) to system size 
and initial energy density has motivated the search for possible alternative non-hydrodynamic 
sources of azimuthal harmonics. 
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Figure 1: Left panel: Experimental v„ harmonics for p+Pb collisions and |A;;| > 2 taken from ATLAS 1121 . The 
magnitudes of the moments are comparable to the one observed in A+A reactions. Right panel: Azimuthal Fourier 
moments arising from the GB distribution (Eq. as a function of the gluon transverse momentum and for dilferent 
choices of the exchanged momentum, as computed in 1221 . 


While it has been shown that hydrodynamic equations, with particular assumptions about 
the initial and freeze-out conditions mi, are sufficient to describe the data, the uniqueness of 
that description is not obvious, especially given the unexpected features of the previously men¬ 
tioned data. One important class of non-hydrodynamic models proposed to explain this puzzle is 
based on the Color Glass Condensate (CGC) and Glasma paradigm involving initial-state non- 
perturbative classical held correlations controlled by a gluon saturation scale, Qsatix,A) IIMIl- 
A simpler perturbative QCD source of multi-gluon azimuthal correlations due to classical non- 
abelian held interference effect was recently presented in GLVB Il22l based on the well-known 
Gunion-Bertsch (GB) LO analytic formula Il2?l . This non-abelian gluon radiation generates non¬ 
zero v„ moments with a shape that closely resembles the experimental one — see Fig. |l(b) The 
magnitudes of those parton level harmonics are however too large and some kind of damping 
mechanism is therefore required. A natural candidate for this job is hadronization. In fact, the 
primary advantage of GLVB, for our purpose of exploring systematic theoretical uncertainties 
associated with hadronization scheme dependence of azimuthal harmonics, is its ease of adapt¬ 
ability to Monte Carlo (MC) multi-particle production in p-i-p, p-i-A and Ah-A collisions via the 
HIJING ll24l type of models as emphasized in ll22]l . 

Thus far CGC/Glasma phenomenology of azimuthal harmonics moments in p-tA and A-hA 
reactions has been limited to the idealized local parton-hadron hadronization scheme ll25l that 
by assumption preserves the parton level correlations into the hnal-state hadrons. This guar¬ 
antees minimal distortion of the initial-state multi-parton correlations predicted by QCD. More 
realistically, it is well-known from decades of hadronic phenomenology that phenomena like 
the production and dacay of intermediate hadron resonances and hnal-state correlations can in¬ 
troduce non-trivial complications and, therefore, some non-peiturbative hadronization scheme 
is required. Hadronization phenomenology cannot be rigorously predicted from QCD, but two 
generic approaches, carefully tuned to e^e^, ep, and pp data, have been developed during the 
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years. These models are the independent fragmentation scheme ll26ll and the LUND string 
model Il27l implemented in the JETSET algorithm Il28l . They can be used to estimate the distor¬ 
tions of initial-state partonic correlations due to a more realistic hadronization process. In par¬ 
ticular, in the nuclear event generator HIJING, JETSET is used to convert the multiple diquark- 
quark beam jets with gluon mini-jet kinks into the multi-hadron resonances with particle data 
book properties and decay branchings. 

In this paper we study the hadronization scheme dependence of final-state correlations using 
three different hadronization models that can be conveniently selected within the JETSET code. 
We deliberately neglect any intermediate parton or hadron level transport or hydrodynamic effect 
to concentrate exclusively on final-state hadronization modification of initial multi-gluon corre¬ 
lations. We compare: (1) local parton-hadron duality (LPH), (2) collinear (independent) parton- 
to-hadron resonance fragmentation/decay (CPR) and (3) LUND string hadronization models. 
We treat Qsat as an input parameter controlling the pj range of azimuthal asymmetry of Gunion- 
Bertsch multi-beam gluon bremsstrahlung. Our results strongly support the analysis of Skokov et 
al. ED, indicating that the CGC/Glasma predictions (without detailed hadronization modelling) 
should be limited to the kinematic range pr > 3 GeV to avoid complications and theoretical 
uncertainty of non-perturbative hadronization physics. 

At this point, it should be stressed that the configurations that we will consider in this work 
are quite simple and hence cannot be considered as a realistic representation of the ridge effect 
observed experimentally in p-i-A reactions. However, as we will see, they are already enough to 
draw strong conclusions about the hadronization mechanism dependence of the final azimuthal 
harmonics. 


2. Elliptic interference harmonics of gluons from two recoiling beam jets 


Eor our simulations, we used the MC hadronization algorithm JETSET 7.4 1^ with 30x10® 
simulated GLVB “events” with two recoiling beam jets. Each beam jet is represented by a high 
invariant mass qq pair along a “beam axis”, t Eor our simulations the invariant mass of each 
beam jet is taken to be 100 GeV. The two beam jets are assumed to scatter with equal but opposite 
net momentum transfer Qi = -Q 2 - {q,if/), with magnitude distributed according to a simple 
Rutherford form ^^d azimuthal direction, ij/, distributed uniformly in 

Eor our numerical simulations, we further assumed Poisson fluctuations of the number, Ng, of 
bremsstrahlung gluons with {Ng) - 8 and distributed uniformly in rapidity between kinematic 
bounds. 

In the LUND color string model Il27l . all gluons are represented by “kinks” of the qq string 
that deform it in the transverse plane. In contrast, from an independent fragmentation point of 
view 1261, gluons are assumed to be just isolated partons that hadronize independently from 
each other. We generate the bremsstrahlung gluons transverse momenta and azimuthal angles, 
k - {k, (p), from the perturbative regularized GB distribution: 




Ql 


q^+p^ 


+ QLf + ^QCd) +q^+H^- 2kq cos(0 - i//)) ’ 


( 1 ) 


with (2sat = 2 GeV being the typical momentum scale expected from the CGC model for p-i-A 
collisions. We take p - 300 MeV and Aqcd = 200 MeV as infrared regulating scales. 

This setup produces initial azimuthal anisotropy at the parton level that with the simplest 
LPH scheme leads to harmonic moments similar to experiment, as seen in Pig{T] In particular, 
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Figure 2: distributions for pairs of initial gluons (blue) and final pions (red) for different pj bins. The hadronization 

has been performed using the independent fragmentation scheme. The solid lines are the fitted curves. 


the pseudo-rapidity independence and the preference of radiated gluons transverse momenta k to 
be aligned along ±Q gives rise to a “ridge-like” structure with azimuthal asymmetry resembling 
the one observed experimentally. 

After the hadronization, an analysis is performed over pairs of final pions with criteria as 
close as possible to ID. In particular, our pseudo-rapidity cuts on the final pions are I 77 I < 2.5 and 
IAt/I > 2 (the so-called long-range), Arj - rja - rji, being the relative psuedo-rapidity of the pair. 
Both pions are taken in the same pj bin and the correlation functions we computed are given by: 


C(A4>; pt) 


S (A 0 ; pt) 
B(A0; Pt) ’ 


( 2 ) 


where A(p - (f>a - (pt is the relative azimuthal angle between the two pions and S and B represent 
the same and mixed event pair distributions respectively 0 . The analysis is repeated for the 
three hadronization schemes to reveal the hadronization scheme dependence of final azimuthal 
harmonics shown in Figs.j^an d 0 

From these distributions we extracted the n < 6 single-particle moments, v^ipr), using the 
following fitting function: 


Cfit(A0; Pt) - Co 


1 +2^vl(pT)cos(nAcf>) 


The results of this fit for n < 4 are shown in Fig. 


(3) 


2.7. Numerical Results 

The computation of the two-particle A(p distributions as reported in Figs. and shows 
that the process of hadronization has some deeply non-trivial consequences for the correlations 
between final state particles. 
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Figure 3: Same quantities as reported in Fig. I^but with hadronization performed using the LUND string model. 


The gluon curves clearly present the behavior expected from the discussion around Eq. Q- 
The long-range near-side (Acp ^ 0) peak is simply due to the fact that the gluons are produced 
with transverse momentum preferentially close to the exchanged momentum Q, while the away- 
side (A(p ^ n) peak is produced by the recoil of the two strings. These features extend over the 
full pt range. 

The distributions of the final pairs of hadrons, instead, have some very peculiar properties, 
depending on which pr region we are considering. For small values of the pion transverse 
momentum — smaller than a scale. A, to be determined in Sec. — the initial anisotropy is 
extremely reduced and the two-pion distributions become more and more uniform for decreasing 
Pt- This behavior is common to both the LUND model and the independent fragmentation 
model, even though the shape of the A0 distributions already appears to be different for pj- > 1 
GeV. The reasons for this strong dilution of the initial anisotropy are proper to the process of 
hadronization itself and will be analyzed again in Sec.|^ If a simple system like ours already has 
such a “decoherence power”, one can reasonably expect the same thing to be valid for a more 
realistic situation — see for example Sec. It should also be stressed that, since the low-px 
pions seem to carry no information about the initial gluon distribution, this feature appears be 
true in general, no matter what causes the partonic correlation in the first place, and hence can be 
applied to other initial-state models as well. 

The second result refers to the higher transverse momentum regime (pj d). In both 
models this region presents strong two-particle correlations among the final pions, since the gluon 
anisotropy is transmitted more efficiently. However, one can immediately appreciate an essential 
difference between the independent fragmentation and the QCD string models. In the former, by 
definition, each parton fragments independently from the others and hence the pion correlation 
function closely resembles the gluon one since, when higher pr are considered, no other relevant 
sources of angular correlation come into play. For the case of the LUND model, instead, while 
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Figure 4: Fitted single-particle azimuthal moments, Vnipr), for initial gluons (upper panel) and final pions as obtained 
from the independent fragmentation scheme (lower left panel) and from the LUND model (lower right panel). The gluons 
do not show any odd Fourier coefficient. This feature is shared by the independent fragmentation case as well while the 
Lund model case presents large negative V] and V 3 . 


the near-side peak is generally lower for the pions than for the gluons, an extremely pronounced 
away-side peak is present. This large Aip ^ n signal (several times bigger than the A0 ^ 0 one) is 
a direct consequence of the transverse momentum conservation taking place at each breaking of 
the color string ll28l . Each qq system, in fact, has zero total pr and the transverse fluctuations are 
governed by a roughly Gaussian probability distribution. This means that when the string breaks 
and a new pair of partons, say q' and is created from the vacuum, the pr of these partons has 
a probability distribution given by 


'P(Pt) exp 


\ 


= exp 



kpI'' 


(4) 


with m being the mass of the produced q' and a: ^ 1 GeV/fm being the string constant. By con¬ 
servation of total momentum they are always produced back-to-back. It then follows that, for 
every hadron of momentum pj there will always be another with momentum -pT, i-e. such that 
A(p - n. This essentially means that the differential distribution C(A0; pj) contains a contribu¬ 
tion proportional to 5(A0 - n), which gets broader when finite pj bins are taken into account, as 
in any realistic situation. We conclude that the QCD strings hadronization scheme that conserves 
event energy-momentum (unlike independent fragmentation) can introduce large away-side cor¬ 
relations in the spectrum of the final hadrons that can lead to negative odd harmonics absent at 
the purely partonic level. 
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Independent fragmentation 

LUND model 

PT [GeV] 

0.25 

0.75 

1.25 

1.75 

0.25 

0.75 

1.25 

1.75 


0.04 

0.10 

0.36 

0.60 

0.06 

0.30 

0.65 

1.00 

KIvl 

0.00 

0.03 

0.25 

0.63 

0.00 

0.20 

0.50 

1.00 


Table 1: Ratios between the pion even moments (v^) and the gluon ones (v,f) for pr <2 GeV. One immediately notices 
that there is a consistent damping of the initial magnitudes. In particular, the final-state are essentially zero (less than 
30% of the initial ones) for pj- < 1 GeV. The suppression for pr > 2 GeV appears to be negligible. 


All these properties lead to the azimuthal moments, v„{pt), shown in Fig. which therefore 
differ dramatically between the two hadronization schemes. In particular, as a consequence of 
the almost complete flatness of the A(f) distributions, both models present even harmonics for 
low-pT which are much smaller than the gluonic ones. To be more quantitative, in Tab. [T] we 
compare the magnitudes of the even azimuthal harmonics at the partonic and hadronic level. As 
one can see, the damping effect is almost total for pr < 1 GeV and persists up to pj- < 2 GeV. 

Moreover, the prominent away-side peak present in the distributions of Fig. causes large 
negative odd harmonics due to assumed local transverse momentum conservation. 

It should also be noted that these v„ arise solely from initial GB gluon bremsstrahlung and 
final hadronization effects. No final state interactions have been included in our simulations. 

3. Hadronization mechanisms that damp gluon harmonics 

At this point one needs to ask; what is the scale which determines the quenching of anisotropy 
for small values of the transverse momenta? Our picture has two natural scales: the typical 
energy of the hadronization process, d ^ 1 GeV, and the typical momentum exchanged, Qsat = 2 
GeV. Which of the two plays the role of discriminant between the uniform and the anisotropic 
regions? 

To answer this question we repeated the previous simulation using the independent fragmen¬ 
tation scheme but with a different value for the average momentum exchanged, Qsat = 4 GeV. 
The results are reported in Fig.|^ If compared to Fig.j^one can appreciate that the final hadronic 
anisotropies now follow the gluonic ones much more closely than in the Qsat = 2 GeV case, as 
one might expect, for almost all values of the transverse momentum except for pr < 1 GeV, 
where the two-pion distributions are again completely flattened out. This shows that the energy 
scale relevant in this decoherence effect is indeed the scale of hadronization, A ^ 1 GeV. 

As a further check, one can study the average ratio between the pion and the gluon momenta 
both for the transverse and longitudinal components. Our simulation indicates that for the two 
models 


(pt) _ [ 0.54 Indep. fragm. {Pz) (2.08 Indep. fragm. 

(^“|o.58 LUND ’ ^“|2.51 LUND ' 

This again is a hint for the role of hadronization in damping the anisotropies at low pr- 

The reasons for this “quenching power” of the hadronization process can be found in es¬ 
sentially two features; (i) non-collinearity of the gluon radiation and (ii) isotropic decays of 
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Figure 5: Same quantities as in Fig.|^and|^but obtained using the independent fragmentation scheme and Qsat = 4 GeV. 
It is evident how increasing the momentum exchanged enhances both the initial and final anisotropies for every value of 
the trasnverse momentum except for pj < A = \ GeV (first two upper panel from the left), where the pions distributions 
are completely uniform. 


resonances. In particular, the fact that gluon fragmentation is not perfectly collinear is already 
manifest in the first order pQCD Altarelli-Parisi distribution. 


dNg 

dzdk^ 


2ti 


Piz), 


(6) 


where k is again the gluon transverse momentum, z is the Bjorken variable and P{z) the splitting 
function. Although this distribution is strongly peaked around k ^ 0 it also has very long tails, 
clearly showing that a totally collinear picture is a too naive approximation (for an interesting 
discussion on the role of fragmentation functions on the azimuthal harmonics see l29l). Mor- 
ever, in the LUND string model we also have the transverse fluctuations of the flux tube, whose 
typical scale is Vic ^ 0.45 GeV (see Eq. Q), which provide another source of dilution of the 
initial anisotropy. Lastly, the intermediate steps between the initial gluons and the final pions 
are populated by resonances. These particles decay without any preferential direction and hence 
strongly contribute to the flattening of the final A0 distributions. In particular, this is the reason 
why the observed decoherence happens for pr GeV, this being the typical mass of the most 
common resonances. To better illustrate this point we performed an overly simplified simulation 
whose description and results are shown in Tab. 

We learn from above that the strong damping of parton level v„ at pj- < T is not a conse¬ 
quence of the particular initial-state model chosen, but it is rather an intrinsic property of the 
hadronization mechanism. This is one of the striking results of our analysis. 
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LUND model 


Independent fragmentation 


Event 1 


Particle 

PT [GeV] 

P, [GeV] 

(j) [rad] 

Particle 

PT [GeV] 

Pz [GeV] 

(p [rad] 

(A++) 

0.90 

0.90 

1.74 

7l~ 

0.22 

0.17 

2.67 

(P°) 

1.00 

0.36 

0.19 

(w) 

1.20 

0.00 

0.09 

(P^) 

0.56 

-0.38 

-0.53 

(P°) 

1.11 

0.91 

-0.02 

P 

0.79 

0.78 

1.95 

tt" 

0.77 

-0.03 

0.01 

7r+ 

0.20 

0.11 

0.79 

71^ 

0.39 

0.11 

-0.74 

7r+ 

0.44 

-0.01 

-0.48 

7r» 

0.11 

0.19 

0.76 

n~ 

0.71 

0.37 

0.58 

71^ 

0.15 

0.19 

0.76 

7r+ 

0.59 

-0.03 

-0.35 

7l~ 

0.29 

0.49 

-0.09 


Event 2 


(P°) 

0.26 

0.08 

0.67 

(P^) 

0.57 

0.36 

0.72 

7 

2.15 

-0.22 

0.04 

(P^) 

0.57 

0.36 

0.72 

7r+ 

0.49 

-0.39 

-1.34 

(P") 

0.81 

0.84 

0.59 


0.42 

-0.24 

-1.55 

(p“) 

0.97 

-0.07 

-0.56 

7r“ 

0.39 

0.10 

3.06 

7r“^ 

0.06 

0.07 

-2.86 

7r+ 

0.33 

0.23 

-0.31 

7T~ 

0.68 

0.35 

-0.92 

71^ 

0.47 

-0.16 

0.54 

7l~^ 

0.19 

0.09 

2.51 

7r~ 

0.21 

0.24 

-2.77 

7r» 

0.64 

0.27 

0.43 

tt" 

0.26 

-0.11 

-2.96 

tt" 

0.44 

-0.07 

1.25 

tt"*" 

0.41 

-0.29 

-3.04 

7l~^ 

1.24 

0.06 

-0.05 


Table 2: Output of two events composed by a qq pair and a gluon for both independent fragmentation and LUND model. 
The quark and antiquark are moving in opposite directions along z-axis, each of them with energy Eq = 10 GeV. The 
gluon, instead, flies away along the jc-axis (?/ = 0, 0 = 0) with energy Eg = 3 GeV. We only report those particles with 
pseudo-rapidity \q\ < 1, i.e. close to the gluon in phase space. One immediately notices that even though the initial gluon 
has zero azimuthal angle, the fragmentation produces particles with (p ^ 0. Every event always contains at least one 
resonance (particles in parenthesis) which then decays isotropically. Moreover, the fact that even the resonances (which 
are produced in the very first step of hadronization) have non-zero anzimuthal angle shows that the gluon fragmentation 
is non-collineai‘. The final pions are therefore widely spread in A0. They also have pr < ^ GeV as expected. 


4. Hadronization damping of triangular and higher gluon harmonics 

As one can see from Eig. in contrast to the LUND model, the independent fragmentation 
leads to very small odd harmonics — albeit fluctuations due to small statistics effects. The 
reason is that for two color antennas the system is back-to-back symmetric, with just two qq 
pairs recoiling from each other. Since at high-pj- no other correlations are introduced, the final 
pions inherit the symmetries of the initial partons. To check the consequences of hadronization on 
the purely geometrical odd moments — and to reproduce a slightly more realistic configuration 
— we performed a simulation involving three color antenna simulating three recoiling beam jets 
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Figure 6: Two-particle distributions obtained from three initial qq pairs in a triangular configuration for both the 
initial-state gluons (blue) and the final pions (red). The hadronization has been performed using the independent frag¬ 
mentation scheme. 


conserving transverse momentum. Both the simulation and the analysis closely follow what we 
explained in Sec.j^but now we implemented a third quark-antiquark pair such that the whole 
system conserves the total momentum, i.e. if 2 1 , 02 and are the momenta exchanged by each 
of the three pairs then Qx+Qi + Qt, -0- In this case we only used the independent fragmentation 
scheme, which lacked the odd v„ in the first place. The results for the A(p distributions and for 
the final Fourier moments are shown in Fig. an d0 

4.1. Results 

As one immediately notices, the initial gluons A(p distributions now have a contribution from 
odd Fourier components, V 2 n+i. Once again, in the pr 4 region none of the initial information 
is preserved and the pions correlation functions are totally flattened. Even though we suffer 
from a lack of statistics, it is also clear that in the higher transverse momentum regime the 
hnal particles distributions keep following the gluonic ones. In terms of Fourier harmonics, 
even though the initial gluons have dehnitely sizeable v„, the hnal pions moments in the low-p/- 
region are essentially zero (less than 10%) — see again Fig.|^ Hadronization is once again very 
effective in quenching the azimuthal asymmetry of the system. For pr ^ 4 the gluonic and 
hadronic v„ have, instead, similar values. 

5. Conclusions 

Our analysis shows that the process of hadronization can lead to major distortions of the az¬ 
imuthal harmonics up to pr - 3 GeV, which are of interest in the search for signatures of perfect 
Huidity in p-i-A systems. In particular, one important aspect seems to be general and model- 
independent: at small values of the transverse momentum of the hadrons the complexity of the 
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Figure 7: Single-particle harmonic moments, v„(pt), for initial-state gluons (left panel) and final-state pions (right 
panels) as extracted from the distributions in Fig. The higher-pj- bins for the pions harmonics clearly suffer low 
statistics fluctuations. 


hadronization itself — namely non-collinearity of fragmentation and isotropic resonance decays 
— greatly reduces the information contained in the initial-state partons, causing an almost total 
smearing of the final hadronic spectrum. As shown in Sec.|^ this is intrinsic to the hadronization 
itself and hence should be taken into account by every model with initial-state anisotropy at least 
with a theoretical error band estimated by testing several hadronization schemes. 

Many models of initial-state anisotropy seem to claim that the only scale relevant for the 
generation of non-zero v„ is the typical momentum exchanged, Qsat- Our analysis shows that this 
is only true above a second hadronization scale, A ^ 1 GeV. Hadronization scheme choice does 
matter, and those schemes that assume local collinear parton-hadron duality or pure collinear 
fragmentation without resonance production may over-estimate final-state hadron harmonics. 
It is interesting to note that assuming that fragmentation is perfectly collinear and described 
by any fragmentation function /(z) should actually lead to an enhancement of the initial gluon 
anisotropy since low-py pions come from gluons with higher transverse momentum that naturally 
have greater azimuthal asymmetry. This is in striking contrast with the results of our JETSET 
simulations, where collinearity of hadronization is broken. 

A second conclusion is that the two-particle correlations in the Pt ^ ^ - 1 GeV region 
can strongly depend on the chosen hadronization model. Specifically, while in the independent 
fragmentation scheme the initial parton anisotropy almost completely transmits anisotropy to the 
final hadrons, in the LUND string model a new source of correlation due to transverse momentum 
conservation is introduced, leading to a large away-side peak in the pion tS.(p distributions that is 
not due to back-to-back mini jets, which are not taken into account in the present simulations. If 
these features also survive to a system more complicated than ours, then this leads to two impor¬ 
tant remarks: first of all, whenever performing a Monte Carlo simulation to study multi-particle 
correlations one should always pay particular attention to the degree of model-dependence of 
the simulation itself since this might introduce important biases on the final conclusions. Sec¬ 
ondly, if one assumes the description of hadronization in terms of a QCD color string as a fairly 
accurate one, then we might have an unexpected source of non-flow away side correlation for 
particles with pr > I GeV which would be difficult to deconvolute from pQCD di-jets by any 
experimental analysis. 

We emphasize again that the properties of the final two-hadron correlations computed in 
this study, for the most part, are due neither to a particular initial-state mechanism nor to any 
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collective transport or hydrodynamic effects; they are genuine consequences of the hadroniza- 
tion process that are scheme dependent. There is unfortunately no guarantee of universality of 
hadronization scheme. 

Lastly, we also notice how, considering the observed mass dependence of final-state correla¬ 
tions Ga, it would be interesting to repeate the previous analysis for different flavors and study 
how hadronization effects are affected by the particle masses. We leave this analysis to a future 
study. 
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